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ABSTRACT 

While  multiresolution  data  analysis,  processing,  and  compression 
hold  considerable  promise  for  sensor  network  applications,  progress 
has  been  confounded  by  two  factors.  First,  typical  sensor  data  are 
irregularly  spaced,  which  is  incompatible  with  standard  wavelet 
techniques.  Second,  the  communication  overhead  of  multiscale 
algorithms  can  become  prohibitive.  In  this  paper,  we  take  a  first 
step  in  addressing  both  shortcomings  by  introducing  two  new  dis¬ 
tributed  multiresolution  transforms.  Our  irregularly  sampled  Haar 
wavelet  pyramid  and  telescoping  Haar  orthonormal  wavelet  ba¬ 
sis  provide  efficient  piecewise-constant  approximations  of  sensor 
data.  We  illustrate  with  examples  from  distributed  data  compres¬ 
sion  and  in-network  wavelet  de-noising. 

1.  INTRODUCTION 

In  the  recent  call-to-arms  found  in  [1],  the  authors  emphasize  that 
spatially  irregular  data  sampling  is  an  inescapable  reality  when 
considering  real-world  sensor  network  deployments.  Using  the  ex¬ 
ample  of  compression  of  a  sensor  network  measurement  field  for 
transmission  to  a  single,  external  sink  via  a  tree-like  routing  struc¬ 
ture,  they  emphasize  that  traditional  regular  2-D  signal  processing 
schemes  simply  do  not  translate  to  the  irregular  setting.  Observing 
that  much  of  sensor  network  signal  processing  literature  assumes 
regularity  of  sensor  placement,  the  authors  motivate  the  need  for 
newer,  irregularity-tolerant  solutions  for  the  sensor  network  set¬ 
ting. 

To  that  end,  we  propose  what  is  to  our  knowledge  the  first 
distributed,  irregular  wavelet  transform  for  sensor  networks.  Em¬ 
ploying  successive  piecewise  constant  approximations  inspired  by 
the  Haar  wavelet  in  the  regular  setting,  our  transform  integrates 
seamlessly  with  existing  hierarchical  routing  schemes  to  align  its 
data  flow  with  efficient  communication  paths  in  the  network.  We 
present  two  variants  on  our  technique  and  outline  exactly  how  to 
implement  each  in  a  real  sensor  network  setting.  Finally,  using  our 
transform,  we  specify  protocols  for  distributed  solutions  to  both 
the  data  compression  and  transmission  problem  discussed  above 
and  the  problem  of  de-noising  sensor  measurements  within  the  net¬ 
work. 

A  good  deal  of  prior  work  has  addressed  wavelet-based  pro¬ 
cessing  in  sensor  networks,  but  none  has  yet  to  resolve  the  diffi¬ 
culties  of  working  with  irregularly-spaced  data  while  appreciating 
the  need  to  minimize  communication  overhead  to  reduce  network 
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power  consumption.  Dimensions  [2,  3]  proposes  an  in-network 
multiscale  wavelet  transform  and  hierarchical  coefficient  routing 
with  its  wavRoute  protocol  but  assumes  sensor  measurements  lie 
on  a  square  grid.  Similarly,  Wisden  [4]  also  assumes  a  square  grid 
and  application  of  regular  wavelets  for  compression  purposes,  as 
do  [5]  and  [6].  Finally,  Fractional  Cascading  [7]  assigns  to  a  sen¬ 
sor  a  view  of  the  entire  network  which  decays  in  resolution  as  the 
network  becomes  more  distant  from  the  sensor.  While  this  ap¬ 
proach  does  not  assume  regular  sensor  placement  and  is  useful  for 
answering  measurement  range  queries,  it  does  not  have  the  broad 
applicability  of  a  wavelet  transform  to  a  variety  of  signal  process¬ 
ing  applications. 

In  specifying  our  transform  algorithms,  we  assume  only  that 
each  sensor  knows  its  location  and  shares  this  information  with  a 
data  sink  and  that  an  efficient  multiscale  routing  hierarchy  is  al¬ 
ready  in  place.  A  number  of  algorithms  for  sensor  localization  ad¬ 
dress  the  former  problem  (see  [8]),  and  both  wavRoute  mentioned 
above  and  COMPASS  [9]  provide  examples  of  the  latter.  COM¬ 
PASS,  in  particular,  creates  a  hierarchical  routing  scheme  by  clus¬ 
tering  nodes,  electing  clusterheads,  and  iteratively  clustering  the 
clusterheads.  Under  this  paradigm,  local  communication  within 
clusters  and  up  and  down  the  routing  hierarchy  are  less  expensive 
than  communication  across  cluster  boundaries.  Adopting  the  rout¬ 
ing  hierarchy  as  our  multiscale  decomposition  hierarchy  allows  us 
to  tailor  data  flows  to  match  the  economics  of  the  routing  protocol. 

In  Section  2  we  review  challenges  of  irregular  data  sampling 
and  motivate  the  need  for  wavelet  transforms  tailored  to  irregular 
data.  Section  3  provides  the  mathematical  foundations  of  the  two 
proposed  wavelet  transforms,  beginning  with  the  multiscale  data¬ 
flow  model  induced  by  the  routing  topology.  The  details  of  the  two 
transforms,  one  a  tight  frame  pyramid  and  the  other  a  complete  or¬ 
thonormal  basis,  follow.  Section  4  addresses  details  of  implement¬ 
ing  these  transforms,  and  Section  5  provides  protocols  for  both  dis¬ 
tributed  data  compression/  transmission  and  distributed  de-noising 
of  sensor  measurements.  Section  6  demonstrates  the  aforemen¬ 
tioned  applications  in  a  simulated  sensor  network  setting,  and  Sec¬ 
tion  7  concludes  with  a  discussion  of  future  directions  for  such 
irregular  wavelet  transforms. 

2.  CHALLENGES  OF  IRREGULAR  DATA  SAMPLING 

A  large  body  of  elegant  wavelet  theory  suited  to  regularly-spaced 
data  has  emerged  in  recent  decades,  but  application  of  wavelet 
analysis  to  irregularly  spaced  samples  is  a  relatively  new  endeavor. 
The  Fourier  mathematics  underlying  regular  wavelets  no  longer 
apply  in  the  regular  setting,  and  second-generation  wavelet  the- 
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Fig.  1.  Voronoi  polygons  (solid)  bounding  a  field  of  sensors  con¬ 
nected  by  edges  in  a  Delaunay  triangulation  (dashed). 


ory  such  as  the  lifting  scheme  [10]  has  risen  to  replace  traditional 
techniques.  Extending  these  ideas  to  two  dimensions  proves  an 
ever  greater  challenge.  The  pyramid-based  approach  of  [1 1]  and 
the  non-redundant  technique  given  in  [12]  represent  some  of  the 
first  sophisticated  attempts  at  tackling  this  problem.  Unfortunately, 
such  approaches  are  typically  intended  to  operate  in  a  centralized 
fashion  on  an  entire  dataset,  so  as  such  they  are  not  directly  appli¬ 
cable  to  the  sensor  network  problem.  For  example,  they  typically 
assume  complete  freedom  to  construct  their  own  multiscale  anal¬ 
ysis  hierarchies,  and  as  stated  in  Section  1,  communication  effi¬ 
ciency  suggests  that  transforms  must  direct  their  flows  along  rout¬ 
ing  hierarchies  formed  with  respect  to  physical  constraints  such  as 
sensor  placement  and  wireless  channel  quality. 

Rather  than  attempting  to  directly  apply  these  ideas,  we  take 
inspiration  from  the  approaches  contained  therein.  Varying  sensor 
density  poses  one  of  the  key  challenges  in  computing  an  irregular 
transform,  as  measurements  from  nodes  in  more  dense  areas  of  the 
graph  will  appear  to  have  greater  “importance”  than  those  from 
less  dense  regions.  To  counter  this  problem,  we  adopt  the  tech¬ 
nique  in  [12]  of  assigning  weight  to  sensor  measurements  based 
on  the  area  of  a  Voronoi  polygon  [13]  around  the  sensor.  Figure  1 
illustrates  a  set  of  Voronoi  tiles  drawn  as  solid  lines  around  a  set  of 
sensor  locations.  The  figure  also  shows  a  Delaunay  triangulation 
of  the  sensor  locations  drawn  in  dashed  lines.  Voronoi  tesselations 
are  in  fact  geometric  duels  of  Delaunay  triangulations,  and  a  sensor 
can  compute  the  area  of  its  surrounding  Voronoi  polygon  merely 
by  knowing  the  edges  for  which  it  is  an  endpoint  in  the  Delaunay 
triangulation.  Delaunay  triangulations  can  be  computed  in  a  dis¬ 
tributed  fashion  within  the  sensor  network,  as  will  be  discussed  in 
Section  4. 

3.  MULTISCALE  TRANSFORMS  WITH  IRREGULAR 
SAMPLES 

We  now  present  the  mathematical  details  behind  our  proposed  trans¬ 
forms.  First,  we  discuss  the  interplay  between  network  routing  and 
multiscale  data  decomposition.  Then,  we  delve  into  the  details  be¬ 
hind  our  two  proposed  transforms,  one  a  tight-frame  Haar  pyramid 
and  the  other  an  orthonormal  basis  Haar  “telescoping”  transform. 

3.1.  Hierarchical  Routing  and  Data  Decomposition 

In  Section  1  we  introduced  the  idea  of  matching  the  transform  hi¬ 
erarchy  to  the  routing  hierarchy.  Such  a  hierarchical  routing  topol¬ 
ogy  will  typically  consist  of  clusters  of  sensors,  one  of  which  is 
designated  as  the  clusterhead.  Figure  2(a)  describes  the  scenario 
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Fig.  2.  Data  flow  (a)  within  a  single  cluster  fci  and  (b)  from  clus- 
terheads  in  fci  through  k$  to  the  supercluster  h  head  in  cluster  kn. 

in  a  typical  6-node  cluster,  labeled  fci.  Under  the  assumed  routing 
economics,  local  communications  within  a  cluster  are  relatively 
cheap,  so  each  sensor  mi  through  m.5  sends  it  measurement  to 
the  clusterhead  rriQ  for  processing.  Using  the  set  m  of  measure¬ 
ments.  the  clusterhead  computes  two  quantities.  A  scaling  coeffi¬ 
cient  Sk-t  describes  some  average  behavior  over  the  cluster,  while 
a  set  dk1  of  wavelet  coefficients  encode  deviations  of  the  measure¬ 
ments  from  the  average  value. 

The  clusterhead  stores  wavelet  coefficients  dk1 ,  but  the  scaling 
coefficient  Sk t  flow  up  to  the  next  level  of  the  transform/routing 
hierarchy,  as  show  in  Figure  2(b).  Along  with  cluster  fci  from  Fig¬ 
ure  2(a),  clusters  through  ki  are  grouped  together  to  form  a 
supercluster  l\ ,  whose  clusterhead  coincides  with  that  of  cluster 
/C4.  Scaling  coefficients  Skx  through  Sk3  are  passed  to  cluster  k,4, 
which  combines  them  with  its  own  Sk4  as  the  input  to  the  next  level 
of  multiresolution  analysis.  A  scaling  coefficient  Six  for  the  super 
cluster  is  generated  along  with  a  set  dix  of  wavelet  coefficients. 
Note  that  the  node  acting  as  the  super  clusterhead  performs  dou¬ 
ble  duty.  As  the  clusterhead  of  cluster  IC4,  it  stores  the  set  dk4  of 
level-fc  wavelet  coefficients,  and  as  the  clusterhead  for  the  super¬ 
cluster  h,  it  maintains  the  set  di4  of  level-1  wavelet  coefficients. 
The  scaling  coefficient  si4  for  the  supercluster  is  then  passed  up 
the  hierarchy  to  participate  in  subsequent  analysis.  The  iterative 
process  terminates  when  one  supercluster  spans  the  entire  network 
and  a  scaling  coefficient  for  a  single  root  node  is  computed.  This 
root  scaling  value,  which  describes  the  average  behavior  over  the 
whole  network,  and  the  entire  set  of  wavelet  coefficients  at  each 
level  represent  the  completed  multiresolution  analysis  of  the  sen¬ 
sor  field.  Note  that,  under  this  model,  the  hierarchical  wavelet 
decomposition  matches  exactly  with  the  routing  topology,  group¬ 
ing  data  arbitrarily  as  dictated  by  routing  concerns.  As  such,  our 
proposed  techniques  are  extraordinarily  flexible  in  the  realm  of  ir¬ 
regular  wavelet  transforms. 


A 


Fig.  3.  Conceptual  mapping  of  piecewise-constant  function  over 
2-D  intervals  of  area  A;  to  piecewise-constant  1-D  function  over 
intervals  of  width  A; 


Given  this  transform  model,  we  can  present  the  details  of  our 
two  proposed  transforms.  Both  are  based  on  piecewise-constant 
approximations  as  seen  in  the  well-known  Haar  wavelet  in  the  reg¬ 
ular  setting  by  the.  We  adopt  a  similar  piecewise-constant  model 
for  the  data  gathered  in  the  network,  assuming  that  measurements 
hold  constant  over  the  Voronoi  polygons  surrounding  each  sen¬ 
sor,  as  described  in  Section  2.  Section  3.2  details  a  pyramid-like 
transform  whose  analysis  and  synthesis  functions  are  remarkably 
easy  to  implement  but  lead  to  a  slightly  redundant  representation, 
computing  1  scaling  and  N  wavelet  coefficients  for  a  cluster  of  N 
sensors.  Section  3.3  presents  a  more  advanced  transform  which 
employs  a  complete  orthonormal  basis  expansion  of  sensor  mea¬ 
surements,  achieving  a  non-redundant  coefficient  representation 
(N  coefficients  for  an  IV-sensor  cluster).  Providing  a  sparser  rep¬ 
resentation  of  the  sensor  field,  this  transform  nonetheless  carries 
a  slightly  higher  cost  in  terms  of  computational  cost  and  ease  of 
implementation. 

3.2.  Tight-Frame  Haar  Pyramid 

The  first  transform,  originally  introduced  in  [14],  represents  an 
extraordinarily  intuitive  and  easy-to-implement  decomposition  of 
sensor  data  into  wavelet  and  scaling  coefficients.  Harkening  back 
to  the  pyramid  image  coders  which  preceded  modern  wavelet  coders 
[15],  the  transform  passes  a  smoothed  average  to  the  next  level  of 
decomposition  via  the  scaling  coefficient  s  and  leaves  behind  as 
many  detail  coefficients  d  as  there  are  sensors  in  a  cluster  —  in 
essence,  a  redundant  transform  with  the  slightest  possible  redun¬ 
dancy. 

We  now  describe  how  to  calculate  the  transform  coefficients 
given  an  TV-length  vector  to  of  measurements  for  a  cluster  of  N 
sensors.  We  assume  that  each  sensor  in  the  cluster  knows  it’s  sup¬ 
port  size  as  defined  by  the  Voronoi  polygon  surrounding  the  sensor. 
These  areas,  {  Ai,  A2, . . . ,  Ajv},  can  be  computed  in  a  distributed 
fashion,  as  detailed  in  Section  4.  Given  the  set  of  areas  and  mea¬ 
surements,  we  can  conceptually  map  the  2-D  transform  problem  to 
1-D  as  shown  in  Figure  3,  where  sensor  measurements  are  taken 
to  be  constant  over  intervals  whose  lengths  correspond  to  sensors’ 
Voronoi  region  areas.  Note  that  this  mapping  is  purely  for  illus¬ 
trative  purposes  and  that  all  subsequent  integrations  in  fact  involve 
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Fig.  4.  Tight  frame  analysis  functions  for  a  cluster  of  5  nodes, 
including  one  scaling  function  S  and  3  wavelet  functions  W\ 
through  W3 


piecewise-constant  functions  over  2-D  regions. 

Given  the  set  {A-,}:^  of  areas,  we  form  the  (TV  +  1)  x  TV 
matrix  K,  given  by 
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where  the  total  area  over  the  cluster  is  given 

as  A tot 

The  bulk  of  the  entries  of  K  are  defined  as 
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while  the  first  sub-diagonal  entries  are  given  by 
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Taken  together,  the  entries  of  K  define  the  basis  functions  given 
in  Figure  4.  The  first  row  of  K  represents  the  constant-valued 
scaling  function,  while  each  ( i ,  j)th  element  in  the  remaining  rows 
gives  the  value  of  wavelet  function  i  —  1  over  the  Voronoi  cell 
surrounding  sensor  j 

Given  the  matrix  K  describing  the  set  of  basis  functions,  all 
that  is  left  to  do  is  simply  multiply  against  the  sensor  measure¬ 
ments  and  integrate  over  the  set  of  support  areas.  Defining  a  an 
TV-by-TV  matrix  A  as  a  diagonal  matrix  of  the  areas  {Aj}!)-!,  we 
form  the  analysis  matrix  as 


Ta  =  KA.  (4) 

Using  the  analysis  matrix,  we  compute  the  ( TV  +  1)  x  1  vector  of 
transform  coefficients  c  as 


(5) 


c  =  T  Arn, 

where  the  first  element  of  c  is  the  scaling  coefficient  s  and  the  re¬ 
maining  N  elements  give  the  wavelet  detail  coefficients  d.  Defin¬ 
ing  the  synthesis  matrix  as 

TS  =  A-'TL  (6) 

we  recover  the  set  of  measurements  from  the  coefficients  as 

m  =  Tsc.  (7) 

As  stated  previously,  the  analysis  and  synthesis  functions  form 
a  tight  frame.  In  fact,  they  form  a  special  class  of  tight  frame 
known  as  a  Parseval  tight  frame  [16],  Given  this  membership, 
we  conclude  that,  even  though  the  frame  is  redundant,  we  have  an 
equivalence  between  energy  in  the  signal  domain  and  energy  in  the 
coefficient  domain: 


N 

E  =  J2  cl  (8) 

i=0 

where  E  is  calculated  by  squaring  the  piecewise-continuous  sig¬ 
nal  given  in  Figure  3  and  integrating  it  over  its  support  regions. 
This  equivalence  implies  that,  were  we  to  selectively  discard  q 
wavelet  coefficients,  the  energy  E  of  the  resultant  signal  would  be 
maximized  by  discarding  the  smallest  q  coefficients  —  exactly  the 
procedure  we  turn  to  for  the  compression  application  discussed  in 
Section  5.1.  Moreover,  given  properties  of  the  transform,  we  can 
bound  the  difference  in  energies  between  the  reconstructed  and 
original  signals  as 

(£-£)  <  2(^c?),  (9) 

i=l  i=l 

ensuring  that  the  approximated  signal  will  not  blow  up  when  the 
inverse  transform  is  applied  —  a  useful  sanity  check  when  working 
with  tight  frames,  whose  redundancy  can  induce  problems  in  such 
situations.  For  a  complete  proof  of  the  transform’s  Parseval  tight 
frame  properties  and  error  bounds,  we  refer  the  reader  to  [17], 
Clearly,  the  tight  frame  pyramid  transform  is  extremely  easy  to 
implement.  Given  a  set  {AjjjLi  of  sensor  support  areas,  comput¬ 
ing  the  wavelet  decomposition  amounts  to  multiplying  the  sensor 
measurements  against  Ta.  which  depends  entirely  on  the  {Aj}^-! 
Unfortunately,  this  does  entail  suffering  redundancy  in  the  coeffi¬ 
cient  representation  —  a  problem  we  can  rectify  with  a  second, 
more  computationally  intensive  piecewise-constant  transform. 

3.3.  Orthonormal  Basis  Haar  Telescope 

Clearly,  a  non-redundant,  distributed  wavelet  transform  for  sen¬ 
sor  networks  is  superior  from  a  data-representation  standpoint.  In 
this  section  we  present  a  technique  which  appeals  to  classic  Haar 
wavelets  in  the  1-D  setting  and  forms  a  set  of  basis  functions  which 
by  construction  are  orthonormal  and  span  the  piecewise-constant 
measurements  in  a  cluster  of  sensors.  Figure  5(a)  recalls  the  reg¬ 
ular  Haar  scaling  and  wavelet  functions  in  1-D.  Applied  to  data 
which  is  piecewise  constant  over  the  intervals  [0, 0.5]  and  [0.5, 1], 
the  scaling  function  merely  sums  measurements  in  the  two  regions 
while  the  wavelet  function  differences  those  measurements.  The 
two  functions  are  each  unit  norm,  mutually  orthogonal  (their  prod¬ 
uct  integrates  to  zero),  and  span  the  entire  interval  (any  function 
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Fig.  5.  (a)  Regular,  1-D  Haar  wavelet  ( W )  and  scaling  (S')  ba¬ 
sis  functions  and  (b)  corresponding  matched  pair  telescope  basis 
functions. 


which  is  piecewise  constant  over  [0,0.5]  and  [0.5, 1]  can  be  formed 
as  a  linear  combination  of  the  Haar  scaling  and  wavelet  functions). 

By  grouping  sensors  in  our  irregular  setting  into  pairs,  we  can 
perform  a  decomposition  mirroring  that  of  the  classic  Haar  case, 
though  we  must  take  into  account  the  support  sizes  of  our  sensor 
measurements  when  doing  so.  Appealing  to  the  1-D  visualization 
of  Section  3,  a  conceptual  mapping  of  sensor  measurements  to  1-D 
intervals  places  paired  sensors  in  adjacent  positions  on  the  line  so 
that  pairwise  basis  functions  such  as  those  shown  in  Figure  5(b) 
apply  over  the  pair  of  intervals.  Again  we  stress  that  this  mapping 
is  merely  for  visualization  purposes  —  the  pertinent  information 
is  contained  in  adjacent  sensor  pairings,  which  we  describe  how  to 
generate  later  in  the  section. 

To  each  pair  of  sensors,  we  apply  the  set  of  scaling  and  wavelet 
coefficients  given  in  the  figure,  where  hi  =  y  Ai(Ai+A2)  ’  = 

-  V  a2(ai1+a2)’  and  =  \J  at W  Thus’  for  a  §iven  cluster’ 

we  have  that 

s  =  (miAi  +  m2A2)y/  A1*A2 

_  _  (10) 

d  =  miAi^/  Al(A12+A2)  +  — m2A2y^  Aa(A1+A2)- 

Using  the  ho,  hi,  and  /i,2  expressions,  it  is  easily  verified  that  the 
pair  of  basis  functions  are  mutually  orthogonal,  have  unit-norm, 
and  span  the  space  of  functions  constant  over  the  Ai  and  A2  in¬ 
tervals. 

Pairing  measurements  to  generate  wavelet  and  scaling  coef¬ 
ficients  according  to  Equation  10  gives  us  roughly  N/2  scaling 
coefficients  over  the  cluster  of  N  sensors.  Recall,  though,  that  to 
fit  into  the  model  of  Figure  2(a)  we  must  describe  the  whole  clus¬ 
ter  with  a  single  scaling  value  for  transport  to  the  next  layer  of 
the  hierarchy.  To  do  so,  we  merely  iterate  the  pairing  and  trans¬ 
form  process,  operating  now  on  scaling  coefficients  whose  sup¬ 
ports  spans  those  of  the  pairs  from  which  they  were  derived.  This 
process  effectively  turns  the  arbitrary,  (N  —  1) — to — 1  single-depth 
tree  provided  by  the  routing  hierarchy  into  a  virtual  binary  tree  as 
shown  in  Figure  6.  We  refer  to  this  structure  as  a  telescope  tree 
since  it  “stretches”  the  cluster  into  a  deeper  tree.  Note  that  the  tree 
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Fig.  6.  Generating  a  virtual  telescoping  tree  from  the  routing  hier¬ 
archy  cluster  specification. 


E-E  =  '£<?, 

i= 1 

where  i  again  indexes  q  smallest-magnitude  wavelet  coefficients, 
which  have  been  set  to  zero  prior  to  reconstruction. 

Computational  complexity  in  this  telescoping  transform  clearly 
resides  in  the  pairing  scheme  used.  We  appeal  to  a  technique 
known  as  Perfect  Matching  [18],  which  pairs  points  so  that  some 
objective  function  is  minimized.  As  clusters  may  have  an  odd 
number  of  nodes,  as  illustrated  in  Figure  6,  care  is  needed  to  spec¬ 
ify  the  non-paired  node,  which  retains  its  value  and  support  into 
the  next  round  of  matching.  We  discuss  the  details  of  matching,  as 
well  as  a  variety  of  other  practical  concerns,  in  Section  4. 


a. 


Fig.  7.  Comparing  a  pair  of  disjoint  wavelet  functions  Wi  and  Wi 
with  the  wavelet  function  W3  and  scaling  function  S  spanning  the 
combined  support  at  the  next  level  of  telescope  hierarchy. 


is,  in  fact,  virtual,  as  the  entire  transform  is  calculated  by  the  clus- 
terhead  given  the  measurements  and  supports  for  all  sensors  in  the 
cluster,  as  described  in  the  introduction  to  Section  3. 

Orthonormality  of  the  basis  functions  holds  over  levels  of  the 
telescope  tree,  as  illustrated  in  Figure  7.  Wi  and  W2  in  the  fig¬ 
ure  describe  wavelet  functions  in  separate  pairs  over  the  intervals 
{Ai,  A2}  and  {A3,  A4}  on  the  same  level  of  the  tree.  Clearly, 
these  functions  are  orthogonal  given  their  disjoint  support  inter¬ 
vals.  W3  and  S  describe  the  wavelet  and  scaling  functions  at  the 
next  level  of  the  hierarchy,  which  pairs  the  scaling  coefficient  over 
{Ai,  A2}  with  that  over  {A3,  A4}.  It  is  not  difficult  to  see  that 
both  W\  and  W2  at  the  lower  level  of  the  hierarchy  are  orthogonal 
to  the  next-level  functions.  By  construction,  all  wavelet  functions 
integrate  to  zero  against  a  constant  value.  As  W\  and  W2  both 
reside  in  regions  spanned  by  constants  in  W3  and  S,  their  inner 
products  with  those  functions  are  zero,  ensuring  orthonormality 
of  the  complete  set  of  wavelet  functions  and  the  scaling  function 
employed  in  the  telescoped  decomposition. 

Given  pairings  that  form  the  binary  tree  in  Figure  6  and  re¬ 
peated  application  of  Equation  10  to  the  pairs,  we  can  form  a  set 
of  N  —  1  wavelets  and  a  scaling  coefficient  for  a  cluster  of  N  sen¬ 
sors.  Moreover,  as  we  have  a  complete  orthonormal  basis  across 
the  cluster,  the  relation  in  Equation  8  and  approximation  energy 
error  is  fixed  at 


4.  IMPLEMENTATION  DETAILS 

Given  both  the  tight-frame  pyramid  algorithm  described  in  Section 
3.2  and  the  orthonormal  basis  algorithm  of  Section  3.3,  we  now 
have  two  methods  for  computing  piecewise-constant  distributed 
wavelet  transforms  for  sensor  networks.  In  this  section  we  discuss 
implementation  issues  of  the  two  transforms,  distinguishing  the 
former  from  the  latter.  But  before  addressing  their  differences,  we 
begin  with  a  common  requirement  of  both. 

As  discussed  in  Section  2,  the  reality  of  irregular  sampling 
requires  that  we  appropriately  weight  sensor  measurements  to  ac¬ 
count  for  non-uniform  sensor  densities.  This  motivates  us  to  calcu¬ 
late  the  area  of  the  Voronoi  polygon  surrounding  each  sensor,  giv¬ 
ing  a  useful  measure  of  the  support-size  of  each  sensor  value.  As 
the  transform  progresses  up  levels  in  the  hierarchy,  support  areas 
are  merely  aggregated  and  applied  to  the  scaling  coefficients  which 
summarize  finer-level  measurements.  Clearly,  all  higher-level  area 
information  can  flow  up  the  hierarchy,  leading  to  efficient  aggre¬ 
gation.  Thus,  the  real  challenge  becomes  efficiently  computing  the 
base-level  Voroni  cell  areas. 

Fortunately,  efficient,  distributed  algorithms  already  exist  for 
computing  a  Delaunay  triangulation  of  sensors.  Recall  from  Sec¬ 
tion  2  that  the  Delaunay  triangulation  is  the  geometric  dual  of  the 
Voronoi  polygon.  Given  the  set  of  Delaunay  edges  to  which  it 
belongs,  a  sensor  can  easily  compute  the  area  of  its  Voronoi  cell 
by  simple  geometric  construction.  A  method  for  computing  the 
triangulation  of  an  n-node  network  with  an  O(nlgn)  commu¬ 
nication  cost  is  given  in  [19],  which  constructs  a  subset  of  the 
Voronoi  graph  having  all  but  the  largest  edges  in  the  triangulation. 
Edges  not  found  by  the  distributed  algorithm  can  conceivably  be 
sent  from  outside  the  network,  as  the  network  user  is  assumed  to 
know  all  sensor  locations  and  can  easily  identify  these  edges.  As 
area  calculation  need  happen  only  once,  this  extra  communication 
overhead  should  have  a  small  amortized  cost. 

Once  we  know  support  sizes,  we  can  implement  each  trans¬ 
form  as  detailed  in  Algorithms  1  and  2,  which  respectively  pro¬ 
vide  pseudocode  for  the  pyramid  and  telescope  transforms  within 
a  cluster.  Note  the  disparity  in  complexity  between  the  two.  The 
pyramid  transform  merely  computes  the  analysis  matrix  specified 
by  Equation  4  given  the  set  of  support  areas  and  multiplies  the 
matrix  against  the  value  set  for  the  cluster,  passing  the  resultant 
scaling  coefficient  and  aggregate  area  up  the  hierarchy. 

The  telescope  transform,  which  yields  a  superior,  non-redundant 
wavelet  expansion,  requires  a  bit  more  computation.  Each  level  of 
the  telescope  tree  within  a  cluster  requires  identification  of  value 
pairs.  To  compute  these  pairings,  we  appeal  to  Perfect  Matching 


Algorithm  1  PyramidXfm( Values,  Areas) 

1:  form  analysis  matrix  T a  (Eqn.  4)  using  Areas 
2:  Coefficients  <—  T AValues 
3:  RETURN  Coefficients!  1),  sum(Areas) 

Algorithm  2  TelescopeXfm! Values,  Areas) 
if  size  of  Values  is  1  then 
2:  RETURN  Values,  Areas 

3:  end  if 

4:  newValues  <—  0,  new  Areas  <—  0 
5:  Pairs  <—  perfectMatch(Areas) 

6:  for  each  pair  €  Pairs  do 

7:  compute  pair  scaling/wavelet  coefficients  (Eqn.  10) 

8:  insert  scaling  coefficient  in  newValues 

9:  insert  sum  of  pair  areas  in  newAreas 

10:  end  for 

11:  if  number  of  values  is  odd  then 
12:  insert  leftover  value  in  newValues 

13:  insert  leftover  area  in  newAreas 

14:  end  if 

15:  TelescopeXfm(newValues.  newAreas) 


as  described  in  [18].  The  Perfect  Matching  algorithm  computes 
the  lowest-cost  match  between  points  in  a  set  given  a  cost  func¬ 
tion.  For  the  telescoping  tree,  we  choose  distance  between  paired 
sensors  as  our  minimizing  metric  and  constrain  chosen  edges  to 
correspond  to  edges  in  the  Delaunay  triangulation  (guaranteeing 
merged  cells  share  a  common  boundary).  As  Voronoi  cells  merge 
into  super-cells,  we  define  the  distance  between  the  super-cells  to 
be  the  minimum  distance  between  centers  of  their  member  Voronoi 
polygons.  This  distance-based  matching  requirement  insures  that 
difference  coefficients  remain  as  small  as  possible  by  exploiting 
the  similarity  in  spatially  proximate  measurements  from  smooth 
fields. 

Though  the  Perfect  Matching  requirement  does  add  computa¬ 
tional  and  implementation  complexity  to  the  transform  computed 
at  a  cluster,  fast  algorithms  exist  for  computing  the  match  [20]. 
Additionally,  the  telescope  tree  within  a  cluster  need  only  be  com¬ 
puted  once  and  stored,  provided  nodes  remain  stationary  and  do 
not  drop  out  of  the  network.  Such  an  approach  also  opens  the 
door  to  computing  all  pairings  outside  the  network  and  flowing 
the  match  information  down  the  hierarchy  to  clusters.  Such  an  ap¬ 
proach  would  entail  far  more  communication,  but  the  amortized 
cost  could  potentially  be  tolerable  when  spread  over  many  subse¬ 
quent  transform  operations. 

Now  that  we  have  completely  described  our  proposed  trans¬ 
forms  and  their  implementation  issues,  we  proceed  in  Section  5 
to  outline  how  transform  coefficients  can  be  applied  to  a  pair  of 
pertinent  signal  processing  problems. 

5.  APPLICATION  OF  TRANSFORM  DATA 

Transmitting  a  compressed  version  of  the  entire  measurement  field 
to  an  exterior  sink  represents  an  expected  sensor  network  task,  as 
suggested  in  [1],  Removing  noise  in  a  distributed  fashion  from 
sensor  measurements  polluted  by  additive  IID  processes  will  also 
find  much  use  in  sensor  networks,  especially  as  a  pre-processing 
step  for  more  advanced  distributed  algorithms.  We  describe  below 
algorithms  for  doing  so,  given  wavelet  coefficients  calculated  as 


Algorithm  3  xmtCompressedfClsthdlD, Threshold) 

1:  CandidateCoeffs  <—  Coefficients(ClsthdlD) 

2:  Coeffs2xmt  <—  0,  ID2xmt  <—  0 

3:  for  each  coefficient  €  CandidateCoeffs  do 

4:  if  (magnitude  of  coefficient)  <  Threshold  then 

5:  insert  coefficient  into  CoeJfs2xmt,  node  ID  into  ID2xmt 

6:  end  if 

7:  end  for 

8:  send  Coeffs2xmt  and  lD2xmt  up  the  hierarchy  to  the  root 
9:  for  each  node  in  cluster  do 
10:  xmtCompressed(  nodelD,  Threshold ) 

11:  end  for 


described  above. 

5.1.  Distributed  Compression  and  Transmission 

Sections  3.2  and  3.3  appeal  to  Equation  8  to  motivate  thresholding 
of  wavelet  coefficients  as  an  optimal  procedure  for  approximating 
the  function  they  represent.  As  reconstruction  error  energy  equals 
that  of  discarded  wavelet  coefficients,  discarding  smallest  coeffi¬ 
cients  first  leads  to  an  optimal  approximation,  given  that  only  a 
specified  number  of  coefficients  can  be  retained. 

Thus,  given  a  magnitude  threshold  below  which  wavelet  coef¬ 
ficients  are  presumed  insignificant,  compression  and  transmission 
of  the  measurement  field  proceeds  as  described  in  Algorithm  3, 
which  is  called  on  the  root  node. 

Given  a  threshold,  a  clusterhead  sends  up  the  hierarchy  sig¬ 
nificant  wavelet  coefficients  along  with  node  identifiers  indicating 
which  nodes  in  the  transform  generated  the  coefficients.  It  then 
passes  the  threshold  down  the  hierarchy  to  child  nodes  in  the  clus¬ 
ter. 

5.2.  Distributed  Measurement  De-Noising 

Noise  in  sensor  measurements  manifests  itself  as  small  wavelet  co¬ 
efficients.  Therefore,  de-noising  proceeds  much  as  compression, 
where  wavelet  coefficients  are  compared  against  a  threshold.  An 
important  difference  arises,  though.  IID  sensor  noise  applies  only 
to  individual  sensor  measurements  and  not  to  their  entire  support 
areas.  As  wavelet  coefficients  are  derived  with  respect  to  continu¬ 
ous  functions  over  the  plane,  we  must  tailor  the  threshold  to  each 
coefficient. 

We  inject  an  estimate  of  of  the  noise  variance  into  the  net¬ 
work.  Considering  the  telescoped  transform  and  a  single  cluster, 
we  fashion  a  threshold  for  each  wavelet  coefficient  in  the  telescope 
tree  as  follows.  Where  {Ai}“=1  and  { Ai}[’_u+1  describe  the  v 
base-level  Voronoi  cells  aggregated  into  the  positive  and  negative 
support  regions  of  the  functions  given  in  Figure  5(b),  the  thresh¬ 
old  t  for  the  wavelet  coefficient  generated  from  those  basis  pairs  is 
computed  as 

U  V 

t  =  On  Ai  +  Aj  (12) 

\  i= 1  j=u+ 1 

Note  that  this  requires  knowledge  at  each  transform  level  of  the 
base-level  Voronoi  cells  spanned  by  the  transform  function  sup¬ 
ports.  Such  information  can  be  computed  in  a  single  bottom-up 
pass  along  the  routing  hierarchy  and  stored  for  future  reference. 
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Fig.  8.  Comparison  of  Haar  Pyramid  and  Haar  Telescope  to  D-2, 
D-4,  D-6,  and  D-8  wavelets  in  a  regular  setting. 


consists  of  a  piecewise-planar  field  with  discontinuities  in  its  left 
and  right  comers,  the  second  of  a  smooth  quadratic,  and  the  third 
of  a  noisy  quadratic  with  a  discontinuity.  Relative  performances  of 
the  Haar  Pyramid  and  Telescope  transforms  are  given  in  Figures  9 
(d),  (e),  and  (f),  with  the  non-redundant  Haar  Telescope  outper¬ 
forming  the  easier-to-implement  Haar  Pyramid,  as  expected.  Fi¬ 
nally  Figures  9  (g),  (h),  and  (i)  show  the  fields  reconstructed  from 
the  Haar  Telescope  function  by  retaining  250  of  the  2,500  coeffi¬ 
cients.  Clearly,  major  features  of  the  fields  are  retained  even  in  the 
presence  of  such  heavy  compression 

6.3.  In-Network  De-noising  of  Measurements 

To  evaluate  the  de-noising  operation,  we  took  samples  from  the 
noisy  quadratic  function  of  Figure  9(c).  Over  60  trials  we  realized 
a  5.48  dB  improvement  in  signal-to-noise  ratio  using  soft  thresh¬ 
olding. 


7.  CONCLUSION 


Now,  given  that  appropriate  thresholds  can  be  applied  to  each 
wavelet  coefficient,  de-noising  proceeds  as  follows.  Starting  at  the 
hierarchy  root,  we  threshold  wavelet  coefficients,  compute  inverse 
transforms,  and  pass  reconstructed  values  down  to  child  clusters 
as  scaling  values  for  subsequent  inverse  transforms.  This  process 
iterates  top-down  until  the  base-level  values  are  recovered.  These 
then  represent  de-noised  versions  of  original  sensor  measurements. 

6.  EXAMPLES 

We  now  proceed  to  demonstrate  the  effectiveness  of  our  proposed 
transforms  on  a  variety  of  simulated  sensor  measurement  fields. 
We  show  results  for  both  compression  and  de-noising,  but  first  we 
provide  a  comparison  of  our  technique,  operating  in  the  regular 
setting,  to  a  standard  regular  wavelet  transform. 

6.1.  Comparison  to  Wavelets  on  Regular  Grids 

As  a  brief  sanity  check,  we  compare  the  performance  of  our  pro¬ 
posed  transforms  in  the  regular  setting  to  transforms  using  regu¬ 
lar  Daubechies-2,  Daubechies-4,  Daubechies-6,  and  Daubechies-8 
wavelets.  For  all  transforms,  we  plot  reconstruction  error  versus 
numbers  of  coefficients  retained  with  increasing  fidelity.  For  this 
setting,  the  sensor  network  consists  of  a  regular  square  grid  of  256 
sensors.  We  allow  sensors  to  cluster  randomly  and  evaluate  the 
resultant  reconstructions  against  those  given  by  the  Daubechies 
wavelets  operating  on  the  entire  square  grid  without  clustering. 
Figure  8  illustrates  the  relative  performances.  Though  both  the 
Haar  Pyramid  and  Haar  Telescope  transforms  are  out-performed 
by  the  Daubechies  wavelets  (as  expected),  their  performance  is 
clearly  comparable.  Thus,  the  transforms  are  capable  of  perform¬ 
ing  well  in  both  the  regular  and  irregular  domains  —  a  claim  the 
Daubechies  wavelets  cannot  make. 

6.2.  Compression  and  Transmission  of  Network  Data 

Now,  we  move  back  into  the  irregular  domain  for  an  examina¬ 
tion  of  the  relative  performances  of  the  Haar  Pyramid  and  Haar 
Telescope  transforms.  Figures  9  (a),  (b),  and  (c)  give  the  three 
simulated  sensor  measurement  fields  used  in  the  study.  The  first 


In  conclusion,  we  have  successfully  demonstrated  novel  irregular 
wavelet  transforms  which  are  ideally  suited  to  the  setting  of  sen¬ 
sor  networks.  Providing  transforms  which  are  both  implementable 
and  practically  valuable,  we  have  demonstrated  their  utility  in  the 
applications  of  distributed  data  compression  and  distributed  de- 
noising.  To  extend  this  work,  we  intend  to  pursue  higher-order 
approximations,  enabling  superior  data  compression. 
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